Lab 7 - More spatial business

Attach packages:

library(tidyverse)
library(tmap)
library(sf)
library(spatstat)
library(maptools)
library(sp)
library(raster)
# library(rgdal)
library(gstat)
library(plotKML) # for points to raster (they won't have this...just testing)

Part 1. Hawaii raster intro

# Read in the raster data

hi_par <- raster("PAR_CLIM_M.tif")
hi_sst <- raster("SST_LTM.tif")
hi_chl <- raster("CHL_LTM.tif")
  
# Base plots
plot(hi_par)

plot(hi_sst)

plot(hi_chl)

First: some useful functions for rasters

Checking it out:

  • crs
  • reprojection
  • cropping
  • simple algebra example
hi_sst@crs # Shows CRS: NAD83
## CRS arguments:
##  +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80
## +towgs84=0,0,0
hi_sst@extent # Shows extent (bounds)
## class       : Extent 
## xmin        : 351124.5 
## xmax        : 959624.5 
## ymin        : 2079231 
## ymax        : 2479731

Example: reprojection to WGS84

wgs84 = "+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs" # Just have this ready to copy/paste

# Reproject
hi_sst_84 = projectRaster(hi_sst, crs = wgs84, method = "bilinear")

# Check the reprojection
hi_sst_84@crs 
## CRS arguments:
##  +proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs +towgs84=0,0,0

raster::aggregate() for resampling

# Sea surface temperature: 
sst_rs <- aggregate(hi_sst, fact = 10)
plot(sst_rs)

# Plot side-by-side for comparison:
par(mfrow = c(1,2))
plot(hi_sst)
plot(sst_rs)

Crop a raster:

# Get these extents from hi_sst_84 (call in console to see) what the actual limits are for hi_sst_84, then decide on cropping boundaries

# First create a spatial polygon
bounds <- as(extent(-156.2, -154.5, 18.7, 20.4), 'SpatialPolygons') # Keep in mind, this could be any polygon shape (state outline, county outline, etc.)

# Reproject
crs(bounds) <- crs(hi_sst_84)

# Then crop: 
sst_crop <- crop(hi_sst_84, bounds)

# And plot:
plot(sst_crop)

A simple algebra example:

Let’s say we’re creating a nonsensical variable called “tropicality”, which is the sum of the PAR + SST + 2*ChlA. How can we create a layer for tropicality?

First let’s reprojeect and get everything into the same CRS:

Use method = “bilinear” for continuous variables, “ngm” for categorical

hi_par_84 <- projectRaster(hi_par, crs = wgs84, method = "bilinear")

hi_chla_84 <- projectRaster(hi_chl, crs = wgs84, method = "bilinear")

# Now we have PAR, Chl-a, and SST all in the same CRS (WGS84) and can start doing some simple algebra.

Plot them side-by-side:

par(mfrow = c(1,3))
plot(hi_sst_84)
plot(hi_chla_84)
plot(hi_par_84)

Raster math is pretty straightforward:

trop <- hi_par_84 + hi_sst_84 + 2*hi_chla_84
## Warning in hi_par_84 + hi_sst_84: Raster objects have different extents.
## Result for their intersection is returned
plot(trop)

We can also explore some stuff about the raster data:

hist(hi_sst_84)

length(hi_sst_84)
## [1] 1020102

And we might want to plot these in tmap instead:

Let’s look at sea surface temperature.

islands <- read_sf(dsn = 'islands', layer = "Island_boundaries") %>% 
  dplyr::select(Island) %>% 
  st_simplify(dTolerance = 10) %>% 
  st_transform(crs = 4326)

plot(islands)

tmap_mode("plot") # or switch to tmap_mode("view")
## tmap mode set to plotting
tm_shape(hi_sst_84) +
  tm_raster(title = "Mean Sea Surface Temperature") +
  tm_layout(bg.color = "navyblue", 
            legend.position = c("left","bottom"),
            legend.text.color = "white", 
            legend.text.size = 0.5) +
  tm_shape(islands) +
  tm_fill("darkgreen")

# Or name it and export
sst_map <- tm_shape(hi_sst_84) +
  tm_raster(title = "Mean Sea Surface Temperature") +
  tm_layout(bg.color = "navyblue", 
            legend.position = c("left","bottom"),
            legend.text.color = "white", 
            legend.text.size = 0.5) +
  tm_shape(islands) +
  tm_fill("darkgreen")

tmap_save(sst_map, "sst.png", height=5)
## Map saved to /Users/ahorst/github/esm-244-lab-7/sst.png
## Resolution: 2251.282 by 1500 pixels
## Size: 7.504274 by 5 inches (300 dpi)

Example: Conditional rasters and masking

Let’s say we have a sensitive species and we’re trying to find suitable habitat. They like warm water (average temp >= 25.6 deg C) and PAR below 54.

# Currently don't have matching extents, we need to update:
extent(hi_sst_84) <- extent(hi_par_84)

# Check compareRaster...nope. Mismatching columns & rows is still a problem. 

# But we also need to make sure they have the same number of rows & columns:
cr <- raster(nrow = 822, 
             ncol = 1229, 
             xmn = -160.4365, 
             xmx = -154.5373, 
             ymn = 18.7309, 
             ymx = 22.44634)

sst_new <- resample(hi_sst_84, cr, method = "bilinear")

compareRaster(sst_new, hi_par_84) # TRUE!
## [1] TRUE

Plot both of them, and crop to a smaller area (for better visualization):

plot(sst_new)

plot(hi_par_84)

Create cropped versions:

# Created 'bounds_main' as earlier: 

bounds_main <- as(extent(-159.9, -159.2, 21.7, 22.3), 'SpatialPolygons') # Keep in mind, this could be any polygon shape (state outline, county outline, etc.)

# Reproject
crs(bounds_main) <- crs(sst_new)

par_kauai <- crop(hi_par_84, bounds_main)
sst_kauai <- crop(sst_new, bounds_main)

# Check out PAR:
plot(par_kauai)

# Then SST:
plot(sst_kauai)

Now we only want to isolate regions where the temperature >= 25.4 and PAR < 54.

# Habitat
par_hab <- par_kauai # just makes a copy
par_hab[par_hab >= 54.0] <- NA

plot(par_hab)

sst_hab <- sst_kauai # also makes a copy
sst_hab[sst_hab < 25.6] <- NA

plot(sst_hab)

par(mfrow = c(1,2))
plot(par_hab)
plot(sst_hab)

So where are the suitable locations where these habitats overlap? raster::mask

suit_hab <- mask(sst_hab, par_hab)
plot(suit_hab)

And make a nice map of the location you’ll recommend:

kauai <- islands %>% 
  filter(Island == "Kauai")

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(suit_hab) +
  tm_raster(legend.show = FALSE) +
  tm_shape(kauai) +
  tm_fill(col = "darkgreen") +
  tm_shape(kauai) +
  tm_borders(col = "yellowgreen", lwd = 2) +
  tm_layout(bg.color = "navyblue")

Part 2. Point pattern analysis

Get the spatial data (counties and red tree voles)

voles <- read_sf(dsn = 'redtreevoledata', layer = "ds033") %>% 
  dplyr::select(COUNTY) %>% 
  filter(COUNTY == "HUM") %>% 
  st_transform(crs = 4326)

# plot(voles)

# Get Humboldt County outline
humboldt <- read_sf(dsn = 'redtreevoledata', layer = "california_county_shape_file") %>% 
  filter(NAME == "Humboldt") %>% 
  dplyr::select(NAME)

st_crs(humboldt) <- 4326

# plot(humboldt)

# Plot them together: 
tm_shape(humboldt) +
  tm_fill() +
  tm_shape(voles) +
  tm_dots(size = 0.2)

# Or with ggplot2: 
ggplot() +
  geom_sf(data = humboldt) +
  geom_sf(data = voles) +
  
ggsave("humvoles.png", 
       units = "in", 
       width = 4, 
       height = 6, 
       dpi = 300)

# Another example (with tiff...there's also jpeg, png, etc.)

# tiff("humvoles2.tiff", units = "in", width = 5, height = 5, res = 300)

ggplot() +
  geom_sf(data = humboldt, fill = "black") +
  geom_sf(data = voles, color = "red", alpha = 0.5)

# dev.off()

We want to explore point patterns in a few different ways. Quadrats. Distance-based methods.

First we need to convert to ‘ppp’ and ‘owin’ - the points and windows, as used by maptools and spatstat (because sf is still catching up for raster and point pattern analysis stuff)

voles_sp <- as(voles,"Spatial")
voles_ppp <- as(voles_sp, "ppp")

humboldt_sp <- as(humboldt, "Spatial")
humboldt_win <- as(humboldt_sp, "owin")

voles_pb <- ppp(voles_ppp$x, voles_ppp$y, window = humboldt_win)
## Warning: 1 point was rejected as lying outside the specified window
## Warning: data contain duplicated points
plot(voles_pb)
## Warning in plot.ppp(voles_pb): 1 illegal points also plotted
vole_qt <- quadrat.test(voles_pb, nx = 5, ny = 10) # nx and ny are number of columns/rows for the rectangles created 
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
# Returns: VoleQT
# Chi-squared test of CSR using quadrat counts

# data:  VolePPP 
# X-squared = 425.94, df = 45, p-value < 2.2e-16
# alternative hypothesis: two.sided 
# Reject the null hypothesis of spatial evenness! But we still don't know if more clustered or more uniform...

plot(voles_pb)
## Warning in plot.ppp(voles_pb): 1 illegal points also plotted
plot(vole_qt, add = TRUE, cex = 0.4)

Plot densities:

point_density <- density(voles_pb, sigma = 0.02)
plot(point_density)

# Can you start viewing this in tmap? Yes, rasterize it: 
wgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
vole_raster <- raster(point_density, crs = wgs84)

# Then plot: 
tm_shape(vole_raster) +
  tm_raster(midpoint = NA, 
            palette = "Blues", 
            legend.show = FALSE)

Nearest neighbor (G-function)

r <- seq(0,0.15, by = 0.005)

gfunction <- envelope(voles_pb, fun = Gest, r = r, nsim = 100, nrank = 2) # Sig level of Monte Carlo = 0.04
## Generating 100 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,  100.
## 
## Done.
plot(gfunction$obs ~ gfunction$r, type = "l", col = "black", lty = 11)
lines(gfunction$hi ~ gfunction$r, type = "l", col = "blue", lty = 8)
lines(gfunction$theo ~ gfunction$r, type = "l", col = "red", lty = 6)
lines(gfunction$lo ~ gfunction$r, type = "l", col = "green", lty = 4)

# Confirms, in combination with quadrat.test, clustered data!

Nearest Neighbor by Ripley’s K (using L standardization)

r2 <- seq(0,0.5, by = 0.05)

lfunction <- envelope(voles_pb, fun = Lest, r = r2, nsim = 20, rank = 2, global = TRUE)
## Generating 20 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,  20.
## 
## Done.
plot(lfunction$obs ~ lfunction$r, type = "l", col = "black", lty = 11)
lines(lfunction$hi ~ lfunction$r, type = "l", col = "blue", lty = 8)
lines(lfunction$theo ~ lfunction$r, type = "l", col = "red", lty = 6)
lines(lfunction$lo ~ lfunction$r, type = "l", col = "green", lty = 4)

Diggle-Cressie-Loosmore-Ford test of CSR

DCLFTest <- dclf.test(voles_pb, nsim = 100, rank = 2) 
## Generating 100 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,  100.
## 
## Done.
DCLFTest
## 
##  Diggle-Cressie-Loosmore-Ford test of CSR
##  Monte Carlo test based on 100 simulations
##  Summary function: K(r)
##  Reference function: theoretical
##  Alternative: two.sided
##  Interval of distance values: [0, 0.250888]
##  Test statistic: Integral of squared absolute deviation
##  Deviation = observed minus theoretical
## 
## data:  voles_pb
## u = 0.001952, rank = 1, p-value = 0.009901

Part 3. Kansas rainfall kriging

# Get Kansas rainfall data
ks_rain <- read_csv("KSRain2.csv")
## Parsed with column specification:
## cols(
##   LAT = col_double(),
##   LON = col_double(),
##   AMT = col_double(),
##   ID = col_character(),
##   ST = col_character()
## )
ks_sf  <-  st_as_sf(ks_rain, coords = c("LON", "LAT"), 
                 crs = 4326)

# Get county data
ks_counties <- read_sf(dsn = 'KSCounties', layer = "ks_counties_shapefile")
st_crs(ks_counties) = 4326

# Plot with tmap:
tm_shape(ks_counties) +
  tm_fill() +
  tm_shape(ks_sf) +
  tm_dots("AMT", size = 0.5)

# Or with ggplot:
ggplot() +
  geom_sf(data = ks_counties, 
          fill = "gray10", 
          color = "gray20") +
  geom_sf(data = ks_sf, aes(color = AMT)) +
  scale_color_gradient(low = "yellow", 
                       high = "red") +
  theme_minimal() +
  coord_sf(datum = NA)

But we want to make predictions across the entire state using kriging.

First, make the rainfall data a Spatial Points data frame:

ks_sp  <- as_Spatial(ks_sf)

Then make a grid that we’ll krige over:

# bbox(ks_sp) to check bounding box of the spatial points
lat <- seq(37, 40, length.out = 200)
long <- seq(-94.6,-102, length.out = 200)

# Then make it into a grid: 
grid <- expand.grid(lon = long, lat = lat)
grid_sf <- st_as_sf(grid, coords = c("lon","lat"), crs = 4326)
grid_sp <- as_Spatial(grid_sf)

Then make a variogram:

# Create the variogram:
ks_vgm <- variogram(AMT ~ 1, ks_sp)

# Look at it: 
plot(ks_vgm)

# Fit the variogram model using reasonable estimates for nugget, sill and range:
ks_vgm_fit <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Sph", range = 200))

# Plot them both together
plot(ks_vgm, ks_vgm_fit) # Cool! So what are the values

# Just FYI: there are other models (Gaussian, Exponential) - how do those line up? 
ks_vgm_gau <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Gau", range = 200))

plot(ks_vgm, ks_vgm_gau)

# You can check the sum of squares of residuals for each: 
attr(ks_vgm_fit, 'SSErr') # 0.00214 (and could compare to other models...)
## [1] 0.002138374
# We'll stick with the Spherical model: 
ks_vgm_fit # Nugget = 0.102, sill = 0.954, range = 235
##   model     psill    range
## 1   Nug 0.1021677   0.0000
## 2   Sph 0.9537363 235.1416

Now, kriging!

ks_krige <- krige(AMT ~ 1, ks_sp, grid_sp, model=ks_vgm_fit)
## [using ordinary kriging]

And visualize it:

ks_krige_df <- as.data.frame(ks_krige) # View it after this to show output

# Rename things to make it a little nicer
ks_krige_2 <- ks_krige_df %>% 
  rename(lon = coords.x1, lat = coords.x2, predicted = var1.pred, err = var1.var)

# Make this into spatial data again
rain_predicted  <-  st_as_sf(ks_krige_2, coords = c("lon", "lat"), 
                 crs = 4326)

# Get Kansas outline to crop: 
ks <- read_sf(dsn = "states", layer = "cb_2017_us_state_20m") %>% 
  dplyr::select(NAME) %>% 
  filter(NAME == "Kansas") %>% 
  st_transform(crs = 4326)

# Crop the rainfall data
rain_cropped <- st_intersection(rain_predicted, ks)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
# Initial plot
plot(rain_cropped) # But this is points

# So is this (cheating...)
# tmap: 
  tm_shape(rain_cropped) +
  tm_dots("predicted", size = 0.05) +
  tm_shape(ks_counties) +
  tm_borders() +
    tm_layout(legend.bg.color = "white", legend.position = c("left","bottom"))